#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<string.h>

#define MAX_PIONT 201

typedef struct 
{
    double x;
    double y;
} Piont;

//生成数据文件
void number_data(const char *some_dat,int t,double a,double b,double c,double d)
{
    FILE *fp=fopen(some_dat,"w");
    if(fp==NULL)
    {
        printf("无法打开文件");
        exit(EXIT_FAILURE);
    }

    for(int i=-100;i<=100;i++)
    {
        double x=i;
        double y=0.0;
        switch(t)
        {
            case 1:y=a*x+b;break;
            case 2:y=a*x*x+b*x+c;break;
            case 3:y=a*x*x*x+b*x*x+c*x+d;break;
            default:
                fclose(fp);
                exit(EXIT_FAILURE);
        }
        y+=(rand()/(double)RAND_MAX)-0.5;
        fprintf(fp,"%lf %lf\n",x,y);
    }
    fclose(fp);
}

//读取数据文件
int read_data(const char *some_dat,Piont* data)
{
    FILE *fp=fopen(some_dat,"r");
    if(fp==NULL)
    {
        printf("无法打开文件");
        return 1;
    }
    int count=0;
    while(fscanf(fp,"%lf %lf",&data[count].x,&data[count].y)==2)
    {
        if(++count>=MAX_PIONT)
            break;
    }
    fclose(fp);
    return count;
}

//最小二乘法
void Least_Squares_Method(Piont *data,int count,int t,double *params)
{
    int n=t+1;
    double A[n][n+1];
    memset(A,0,sizeof(A));
    for(int i=0;i<n;i++){
        for(int j=0;j<n;j++){
            for(int k=0;k<count;k++){
                A[i][j]+=pow(data[k].x,i+j);
            }
        }
        for(int k=0;k<count;k++){
           A[i][n]+=data[k].y*pow(data[k].x,i); 
        }
    }
    for(int k=0;k<n;k++){
        int max_row=k;
        //高斯消元
        for(int i=k+1;i<n;i++){
            if(fabs(A[i][k])>fabs(A[max_row][k])){
                max_row=i;
            }
        }
        for(int j=k;j<=n;j++){
            double temp=A[k][j];
            A[k][j]=A[max_row][j];
            A[max_row][j]=temp;
        }
        for(int i=k+1;i<n;i++){
            double factor=A[i][k]/A[k][k];
            for(int j=k;j<=n;j++){
                A[i][j]-=factor *A[k][j];
            }
        }
    }
    for(int i=n-1;i>=0;i--){
        params[i]=A[i][n];
        for(int j=i+1;j<n;j++){
            params[i]-=A[i][j]*params[j];
        }
        params[i]/=A[i][i];
    }
}
//剩余平方和
double squares_sum(Piont *data,int count,double *params,int t)
{
    double R=0.0;
    for(int i=0;i<count;i++){
        double predication=0.0;
        for(int j=0;j<=t;j++){
            predication+=params[j]*pow(data[i].x,j);
        }
        R+=pow(predication-data[i].y,2);
    }
    return R;
}

int main(int argc, char* argv[]){
    if (argc < 2) {
        printf("用法:\n");
        printf("生成数据: %s gendat filename t a b [c [d]]\n", argv[0]);
        printf("曲线拟合: %s fitting filename t\n", argv[0]);
        return 0;
    }
    
    if (strcmp(argv[1], "gendat") == 0) {
        if (argc < 6) {
            printf("参数不足\n");
            return 0;
        }
        
        char* filename = argv[2];
        int t = atoi(argv[3]);
        double a = atof(argv[4]);
        double b = atof(argv[5]);
        double c = (argc > 6) ? atof(argv[6]) : 0.0;
        double d = (argc > 7) ? atof(argv[7]) : 0.0;
        
        number_data(filename, t, a, b, c, d);
        printf("数据已生成并保存到 %s\n", filename);
        
    } else if (strcmp(argv[1], "fitting") == 0) {
        if (argc < 4) {
            fprintf(stderr, "参数不足\n");
            return 1;
        }
        
        char* filename = argv[2];
        int t = atoi(argv[3]);
        
        Piont data[MAX_PIONT];
        int count = read_data(filename, data);
        
        double params[4] = {0}; // 最多支持三次函数
        Least_Squares_Method(data, count, t, params);
        
        printf("最小二乘法拟合结果:\n");
        for (int i = 0; i <= t; i++) {
            printf("参数%d: %.6f\n", i, params[i]); // 提高显示精度
        }
        
        double R = squares_sum(data, count, params, t);
        printf("剩余平方和(R): %.6f\n", R);
        
    } else {
        printf("错误 \n");
        return 1;
    }
    
    return 0;
}